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ABSTRACT 

We present results from ab initio simulations of liquid water-hydrogen mixtures in the range 
from 2 to 70 GPa and from 1000 to 6000 K, covering conditions in the interiors of ice giant planets 
and parts of the outer envelope of gas giant planets. In addition to computing the pressure and 
the internal energy, we derive the Gibbs free energy by performing a thermodynamic integration. 
For all conditions under consideration, our simulations predict hydrogen and water to mix in all 
proportions. The thermodynamic behavior of the mixture can be well described with an ideal 
mixing approximation. We suggest a substantial fraction of water and hydrogen in giant planets 
may occur in homogeneously mixed form rather than in separate layers. The extend of mixing 
depends on the planet’s interior dynamics and its conditions of formation, in particular on how 
much hydrogen was present when icy planetesimals were delivered. Based on our results, we do 
not predict water-hydrogen mixtures to phase separate during any stage of the evolution of giant 
planets. We also show that the hydrogen content of an exoplanet is much higher if the mixed 
interior is assumed. 

Subject headings: Physical Data and Processes: equation of state; planets and satellites: gaseous planets; 
planets and satellites: Jupiter, Saturn, Uranus, Neptune. 


1. Introduction 


The past decade has been a period of extraor¬ 


dinary discoveries in the field of exoplanets (Ex¬ 


oplanet Arch. 20141. For the first time, we now 


have quantitative estimates for occurrence rates of 
different types of planets in our galaxy. A recent 


study by Petigura et al. (20131 that focused on 


planets with periods up to 100 days in the Kepler 
sample demonstrates the prevalence of Neptune 
and sub-Neptune exoplanets. In this sample, the 
planets with a radius larger than 1.5 Earth radii 


have a mean density (Weiss & Marcy 20141 that 


implies that they are composed of both heavy el¬ 
ements (rocks and ices) as well as gas (hydrogen 


and helium) (for mass-radius relation, see Seager 


et al. (2007)). Despite the absence of sub-Neptune 


planets in our solar system, we may be able to 


place constraints on the formation process of ice 
giants by studying the interior and the evolution 
of Uranus and Neptune. Both planets have sizable 
gaseous envelopes surrounding cores composed of 
heavier elements. 


The core accretion model ([Pollack et al. 1996 


Helled &: Bodenheimer|[T014[ [Bodenheimer fc Lis-j 


sauer 2014) suggests that the giant planets form 


in two distinct phases. First rocky and icy plan¬ 
etesimals accumulate to form a dense core. Once 
a critical core mass has been reached, a runaway 
accretion of hydrogen-helium gas sets in and lasts 
until all gas from the planet’s neighborhood has 
been depleted. Most often it has been assumed 
that very little gas is present when the initial core 
forms (Pollack et al. 1996). Because of this as¬ 


sumption, planetary interior models typically as¬ 
sume a dense core of rock and ice with a sharp 
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transition to a gaseous outer envelope that is com¬ 
posed of hydrogen, helium, and a small fraction of 


heavier elements (Stevenson|1985| Hubbard et al. 

1995 Hubbard 119991 

Guiliot|19991 

20051 Militzer 

et aZ.|2008 

). While the heavier element fraction in 


Jupiter’s envelope has been measured in situ by 
the Galileo entry probe (Wong et al.||2004), there 


is no direct measurement that characterizes the 
state of a giant planet core. It is thus plausible 
that some mixing of gas and icy planetesimals has 
occurred when a giant planet formed and that it 
may have persisted until today. 

The mean density of Uranus and Neptune, on 
the other hand, suggests that they have a much 
higher fraction of heavy elements than Jupiter 


and Saturn (Hubbard 1984 pp. 282-295). Based 


on cosmological abundances, it is assumed that 
water is among the dominant species, followed 
by methane and amonia, even though only trace 
amounts of it have been detected spectroscopically 


naz 


in the atmospheres of Uranus and Neptune (Encre- 
]|2003). This paper focuses on water-hydrogen 


mixtures because of the prevalence of water in our 
solar system and the possibility that gases and ices 
mix when giant planets form. 

A recent study (Fortney et al. focused on 

the late accretion of planetesimals when a large, 
gaseous envelope is present. It was predicted 
that even large planetesimals of 100 km in diam¬ 
eter may disintegrate upon entry and the mate¬ 
rial would be distributed throughout the gas en¬ 
velope. Depending on the mixing properties of 
hydrogen and heavier elements, this material may 
either remain in the envelope or gradually settle 
onto the existing core. We will study hydrogen- 
water mixture in this article because little infor¬ 
mation is available about how well hydrogen mixes 
with other elements at high pressure. The phase 
separation of hydrogen and helium has been pre¬ 


dieted to occur in the interiors of Saturn (Steven- 

son & Salpeter 

1977 

Fortney & Hubbard 2004[ 

Morales et al.|2009 

Soubiran et al.|2013) and also 

of Jupiter (Wilson & Militzer|2010 

). The resulting 


release of gravitational energy has been named as 
a reason for the observed excess in Saturn’s lumi¬ 
nosity. It would be of interest to know whether a 
similar process involving the separation of water- 
hydrogen mixtures could operate in the interiors of 
ice giant planets. If a phase separation occurs, the 
envelope would be progressively depleted of heavy 


elements, which affects the interior structure and 
the luminosity of a giant planet. 


On the contrary, it is possible that the core 
of an initially differentiated planet would dissolve 
into the surrounding layer of hot, dense hydro¬ 
gen. Results from recent ab initio simulations 
predict that the rocky and icy components of the 
cores in Jupiter and Saturn are miscible in metal¬ 


lic hydrogen (Wilson & Militzer 


al. 2013 Gonzalez-Cataldo et al. 


2012a Wahl et 


20141. While a 


mixture state is thermodynamically preferred, it is 
not known how fast the cores of giant planets erode 
because gravitational forces counteract the advec- 
tion of heavy elements. This may give rise to a 
semi-convective regime that has already been sug¬ 


gested to occur in the gas giant interiors (Leconte 
fc Chabrier|[M^[M^ . 


Using recent updates of the gravitational mo¬ 


ments of Uranus and Neptune (Jacobson 2007 


20091, revised interior models have been con¬ 


structed. Helled et al. (2011) fitted a density 


profile to these data and showed that they could 
be satisfied by a single layer of hydrogen and he¬ 
lium with a compositional gradient of heavier el¬ 
ements. On the other hand, [Nettelmann et ^ 
(2013) matched the available constraints by con¬ 
structing a typical three layers model assuming a 
rocky core, surrounded by an intermediate layer 


of water in a liquid or superionic state (Cavaz- 


zoni et ^|1999[ [French et a/.||200'9al [Wilson ^ 


al. 


2013), and a gaseous outer envelope. While 


Helled’s model assumes water and hydrogen are 
completely miscible throughout the planet inte¬ 
rior, Nettelmann et al. conversely assume both 
fluids would not mix at the interface of the two 
layers at approximately 10 GPa and 2000 K. This 
underlines why the mixing properties of water and 
hydrogen are important for giant planet interiors. 


Both aforementioned models assume an adi¬ 
abatic, fully convective behavior for each layer. 
This assumption is not consistent with the ob- 


served luminosity of Uranus ( 

Pearl et al. 

1990 

Pearl & Conrath|1991 ) 

Podolak et al. 

(1990 

) sug- 


gested a stably stratified interior model instead. 
Nevertheless, a convective layer of a conducting 
material is needed to sustain a magnetic dynamo 
and to produce the quadrupolar fields observed for 


Uranus and Neptune (1 

^ess et a/.|1986 

1989 

Gon- 

nerney et al.[[1987 

199 

L). Water is assumed to be 


the dominant species in the layer where magnetic 
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fields of ice giant planets are generated. If wa¬ 
ter and hydrogen mix, a gradient of composition 
may be introduced into the planet’s interior dur¬ 
ing formation. It would thus be possible for semi¬ 
convection to operate in ice giant planets which 
would reduce the overall heat flux and allow for a 
vigorous convection in some layers. 


Some experimental data for water-hydrogen 
mixtures are available. Seward & Franck (1981) 
studied mixtures up to 0.25 GPa and observed 
a phase separation below 650 K. Based on this 
result, the presence of water clouds in the deep 
atmosphere of the ice giants has been inferred 


Fegley & Prinn 

1986 

). More recently. 

Bali et 


al. (2013) investigated the properties of water- 


hydrogen in the presence of minerals under deep- 
Earth conditions up to 2.5 GPa. Hydrogen was 
released due to chemical reactions between water 
and the minerals. The analysis of micro-inclusions 
of fluid within the minerals showed that phase 
separation of hydrogen and water was possible at 
temperatures below 1200 K. These findings fa¬ 
vor differentiated interior models for Uranus and 
Neptune. 


At the present time, no experimental data 
are available for higher pressure and tempera¬ 
ture but materials under such conditions can be 
studied efficiently with the ab initio simulations 
that we used throughout this paper. Results from 
such simulations have been shown to agree very 
well with shock wave measurements for hydro¬ 


gen ( 

Lenosky et al. 

1997 

Knudson et al. 

2004 

Loubeyre et al. 

2012 

I. Similarly results of shock 


wave experiments of water have been matched 
with ab initio simulations (Knudson et al. 2012). 
Recently improvements have been made to com¬ 
pute the dielectric constant of water with such 
methods ( Pan et ^|2013 ). Here we use the same 
technique to study water-hydrogen mixtures from 
2 GPa to 70 GPa and from 1000 K to 6000 K. For 
various mixing ratios, we compute the equation 
of state, in particular the pressure and internal 
energy as function of density and temperature. In 
addition we derive the Gibbs free energy of mixing 
by performing a thermodynamic integration. We 
show that the mixtures behave close to an ideal 
mixture. The entropy of mixing is the dominant 
term in the Gibbs free energy of mixing indicating 
that no phase separation occurs at all conditions 
under consideration. 



Fig. 1.— Pressure-temperature diagram with the 
predicted interior profiles for the solar giant plan¬ 
ets. The different expected phase transitions are 
also plotted. The light brown region is where Bali 
et al. found a phase separated system. The cyan 
region shows the parameter range studied in this 
work. There, no phase separation was found. The 
plus symbols mark specihc simulation conditions. 


We also analyze changes in ionic species that 
are present in the mixture as a function of the pres¬ 
sure and the temperature. Finally we comment on 
some implications of computed miscibility proper¬ 
ties. We show for instance that a sharp transi¬ 
tion from a hydrogen-rich envelope to a water-rich 
phase is thermodynamically unstable and that a 
three layer picture for the icy giants is a simplifica¬ 
tion. Furthermore we compare the mass-radius re¬ 
lationship for Neptune-like and sub-Neptune plan¬ 
ets depending on their differentiation. We found 
that for given radius and mass, a fully differen¬ 
tiated planet has a much lower hydrogen content 
than a homogeneously mixed planet. 


2. Simulation methods 


Our investigations of the hydrogen-water mix¬ 
ture rely on ab initio molecular dynamics (MD) 
simulations in which the nuclei are treated as 
classical particles while the electrons are consid¬ 
ered quantum mechanically using density func¬ 


tional theory (DFT) (Hohenberg & Kohn 1964). 
We used the Vienna Ab initio Simulation Pack¬ 


age (VASP) (Kresse & Furthmiiller 1996). The 
time step of the MD simulations was set to 0.2 fs, 
which is short enough to describe the molecular 
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vibrations accurately. Each simulation was per¬ 
formed for a minimum of 0.5 ps and up to 7 ps 
for the lowest temperatures. To keep the temper¬ 
ature constant, we employed a Nose thermostat 
(Nose 1984 1991). For the electrons, we used a 


Fermi-Dirac distribution within a finite tempera¬ 
ture scheme ( Mermin| 1965). We used projector 
augmented wave (PAW) pseudopotentials (Blochl 


1994) with a cut-off radius of = 0.8 Oq for 
hydrogen and = 1.1 oq for oxygen. We used 
a plane-wave basis energy cut-off at 1100 eV. To 
sample the Brillouin zone, we used the Baldereschi 
point (Baldereschi 1973), except for pure water. 


for which a 2 x 2 x 2 Monkhorst-Pack K-point 
grid ( Monkhorst fc Pa^|1976 ) was used. We per¬ 
formed a few simulations with finer K-point grids 
and found consistent results within the statisti¬ 
cal errorbars. For the exchange-correlation func¬ 
tional, we chose the generalized gradient approx¬ 
imation (GGA) of the Perdew, Burke & Ernzer- 


hof (PBE) type (Perdew et al. 1996) because it 


gave reasonable results for pure hydrogen (Gaill- 


abet et al. 2011 Loubeyre et al. 

2012 

) and pure 

water ( 

French & Redmer 

2009b 

Knudson et al. 


2012). We also verified our predictions by per¬ 


forming additional simulations with the van der 


Waals density functional (vdw-DF) by Dion et al. 
(2004); Klimes et al. (2011). The resulting mixing 


properties were in agreement with our PBE pre¬ 
dictions at 1000 K. A more detailed comparison of 
the PBE and van der Waals functionals is given in 


the recent work by Santra et al. (2013). 


To explore the phase separation at a given pres¬ 
sure and temperature, it is necessary to determine 
the Gibbs free energy. Standard molecular dynam¬ 
ics only provide the internal energy and the pres¬ 
sure but not the entropy of the system. Therefore 
we performed a thermodynamic integration using 


an auxiliary classical pair potential {|de Wijs et al. 

1998||Morales et a/.|2009||Wilson & Militzer 

2010 

2012a|b| McMahon et al. 

2012 

Wahl et al. 

2013) 


because it provides the Helmholtz free energy dif¬ 
ference Fi — Fq between two systems characterized 
by two different potentials C/o({r}) and {7i({r}): 


Fi-Fo= [ {Ui-Uo)xdX, (1) 

Jo 

where the parameter A defines a hybrid poten¬ 
tial U\ = Uo + X{Ui+Uo). The {.)\ means that we 
take the average over a trajectory computed with 


the potential U\. We performed the integration 
in two distinct steps. First we integrated between 
the potential as given by the DFT, GdfT) and a 
set of classical pair potentials, Ud, that we con¬ 
structed by matching the forces on configurations 
taken from a DFT trajectory that was computed 
beforehand for a particular temperature, density, 
and composition. 


While pair potentials provide a sufficiently 
good description of water at high pressure for 
our thermodynamic integration technique to work 


well (Wilson & Militzer 2012a Wilson et al. 2013), 


special care needs to be taken to describe molec¬ 
ular hydrogen. Pair potentials have a deep mini¬ 
mum that represents the intramolecular binding. 
If such attractive pair potentials are used with¬ 
out any repulsive many-body term, unphysical 
chains and clusters of hydrogen nuclei form. In 
Militzer (2013), a repulsive many-body potential 


was constructed and a stable thermodynamic inte¬ 
gration technique was obtained for molecular and 
partially dissociated hydrogen. Since here we are 
dealing with water-hydrogen mixtures that lead to 
a more diverse set of short-lived chemical species, 
we pursued a different and simpler approach. We 
constructed sets of non-honding pair potentials 


that we have developed and tested recently (Wahl 
et al. ]M^ . By removing all attractive parts from 
the pair potentials, we prevented the formation of 
unphysical clusters. Hydrogen molecules still form 
gradually as we switch from the non-bonding po¬ 
tentials to the DFT forces. We determined that 
using between 5 to 7 A points was still sufficient to 
accurately calculate the integral in eq. ([^. Exam¬ 
ples of the evolution of this average as a function 
of A are shown in Fig. 


To derive the Helmholtz free energy of the clas¬ 
sical system, Fd, we performed another thermo¬ 
dynamic integration to a reference system with 
known Helmholtz free energy, Fq. For this pa¬ 
per, we used the ideal gas as reference. For this 
integration, we used many more A-steps as it only 
requires classical simulations, which are 10® times 
faster than DFT computations. 


With this two-steps integration approach, we 
are able to derive the Helmholtz free energy for 
the DFT system at different densities, tempera¬ 
tures, and concentrations. Adding the PV term, 
we obtain the Gibbs free energy, Gdft = Adft + 
PdftI^j where Pdft is the pressure given by the 
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DFT calculation and V the volume of the simula¬ 
tion cell. 




Fig. 2.— (C/dft —t^ci) vs A for three different P-T 
conditions for a mixture ratio iVH 2 o:-^H 2 =24:48. 

We considered different concentrations of water 
molecules, x, 


-A^HaO 


( 2 ) 


where W is the number of molecules of type i in 
the simulation cell. We performed simulations for 
the following ratio: 48:0, 40:16, 32:32, 

24:48, 16:64, and 0:92. Except when noted oth¬ 
erwise, the concentrations refers to the total con¬ 
tents in the simulation cell regardless what chem¬ 
ical species form at various pressures and temper¬ 
atures. In Fig. we show a snapshot from a sim¬ 
ulation for X = 0.5 and A^HaO^-A^Ha =32:32, which 
implies that 32 oxygen and 128 hydrogen atoms 
were present in the simulation cell. When we var¬ 
ied the water concentration, we always replaced 
one water molecule with two hydrogen molecules 
so that the volume of the simulation cell would 
not change too much. We performed simulations 
in pressure range from 2 to 70 GPa along four 
isotherms at 1000, 1500, 2000 and 6000 K. 


3. Results 


Fig. 3.— Snapshot of a simulation at 6000 K 
and 45 GPa for a A^H20'-A^H2=32:32 mixture. The 
isosurfaces show the electronic density. The bond 
structure is based on the nuclei distances. 



Fig. 4.— Density versus pressure at 6000 K for 
different water concentrations, x, indicated in the 
legend. The squares represent our simulation re¬ 
sults. The full lines are spline interpolations for 
the pure systems while the dashed lines are the 
predictions from an ideal mixing approximation 
that very well reproduces our direct simulation re¬ 
sults of mixtures. 


From the DFT molecular dynamics and the 
thermodynamic integration, we extracted pres¬ 
sure, internal energy, and Gibbs free energy as 
function of temperature, density, and concentra¬ 


tion. The results are given in Tab. ][^ 

In Fig. |4]|^ we plotted various thermodynamic 


^The table is available in the published article in ApJ. 
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Fig. 5.— Internal energy versus pressure for dif¬ 
ferent concentrations x reported in Fig. The 
dashed lines are again the predictions from an 
ideal mixing approximation. 



Fig. 6.— Gibbs free energy versus pressure for 
different concentrations x reported in Fig. 


quantities for different pressures and concentra¬ 
tions at 6000 K . We used a spline interpolation 
for the pure water and pure hydrogen system. We 
compared our simulation results for the mixtures 
with an ideal mixing approximation using our sim¬ 
ulation results for H 2 and H 2 O and an additive 
volume law at constant pressure and temperature. 
The ideal mixing approximation was found to re¬ 
produce our simulation results for the mixtures 
well. Some small deviations of the order of a few 
percents (up to 10% locally), in particular for the 
density and for the internal energy, can be identi¬ 
fied, however. Nevertheless we conclude that the 


ideal mixing approximation is robust and will be 
sufficiently accurate for the construction of most 
planetary interior models. 

We computed the Gibbs free energy of mixing: 

AG(a^,P,T) = G{x,P,T)-xGn,o{P,T) 

- (l-cc) GH,(P,r), (3) 


where G is a Gibbs free energy per molecule, and 
Gi is the Gibbs free energy of the pure system of 
molecule i. A homogeneous mixture is the ther¬ 
modynamically preferred state if: 


a^AG 

dx"^ 


> 0 , 

P,T 


(4) 


which implies that AG is a convex function of x 
at given pressure and temperature. If the Gibbs 
free energy difference shows a partial or a full con¬ 
cavity, the homogeneous system is unstable and 
undergoes a partial or complete phase separation. 

In Fig. we plotted the Gibbs free energy of 
mixing for different temperature and pressure con¬ 
ditions. Most of the curves are convex and when 
they are not, a convex curve can still be drawn 
within the errorbars. We therefore conclude that, 
for the whole set of conditions we explored, a ho¬ 
mogeneous mixture of hydrogen and water in any 
proportion is the thermodynamically stable state. 
We were not able to identify an indication for 
a phase separation at any condition that we ex¬ 
plored. 

For given pressure and temperature, we can 
split the Gibbs free energy of mixing into three 
terms, AG = AE + PAV — TAS. The inter¬ 
nal energy term, AE, represents the interaction 
between the different species in the mixture. The 
PAV term measures deviations in density from an 
ideal mixture. Finally, —TAS is the full entropy 
of mixing including ideal and non-ideal contribu¬ 
tions. 

From the example at 6000 K and 70 GPa shown 
in Fig. we can infer that the contributions to 
the Gibbs free energy of mixing from the inter¬ 
nal energy and the PAV term are small while the 
entropy is the dominant term by far. The com¬ 
puted entropy can be very well approximated by 
the ideal entropy of mixing of water and hydrogen 
molecules Sid = —/cb [a^lnx-I-(1 — a:) ln(l — x)]. 
We find this approximation to work very well even 
if the system is almost fully dissociated at, e.g., 
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Fig. 7.— Gibbs free energy of mixing as a func¬ 
tion of the water concentration x (hull diagram) 
at different temperatures and pressures. 


6000 K and 70 GPa (see Fig. |^&[^. If one tries 
to use the entropy of mixing for two atomic sys¬ 
tems instead, the agreement with the simulation 
results is inferior (Fig. [^. We conclude that at 
high temperature when many short-lived species 
are present, the attraction between ionic species is 
still sufficiently large so that molecular entropy of 
mixing is a much better approximation than rely¬ 
ing on a mixture of atoms. 


We performed a systematic study of the dif¬ 
ferent contributions in the Gibbs free energy of 
mixing and always found that the entropy is the 
dominant term that can be matched well with an 
ideal mixing approximation for molecules. The 
largest deviations from this approximation arise 
for high water concentration. The presence of 
hydrogen appear to slightly alter the dissocia¬ 
tion fraction of water molecules in the mixtures. 
This is the strongest non-ideal mixing effect that 
we identified. Similarly, the presence of helium 
atoms appears to increase the stability of hydro¬ 
gen molecules when various hydrogen-helium mix¬ 
tures are compared for given pressure and temper¬ 
ature ( Vorberger et a/.|[2007 ). 



Fig. 8.— Decomposition of the Gibbs free energy 
of mixing as a function of the water concentration 
X at 6000 K and 70 GPa. The black dotted line 
shows the predicted contribution for an ideal en¬ 
tropy of mixing of hydrogen and water molecules 
while the dash-dotted line is an ideal entropy of 
mixing of oxygen and hydrogen atoms. The latter 
does not agree well with our simulation results. 


In addition to studying various thermodynamic 
functions, we also determined the different chem¬ 
ical species present in the mixture for the dif¬ 
ferent pressure, temperature, and concentration. 
We used a similar species analysis method as was 
employed by Vorberger et al. (2007). At each 
time step, we computed the distances between 
the nuclei. If a pair of nuclei remained closer 
than a given distance for a minimum duration (ten 
times the vibration period of molecular hydrogen 
th 2=7.6 fs) we considered the two nuclei bound. 
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Fig. 9.— Distribution of the oxygen atoms 
among the different oxygen bearing species at 
2000 K (□) and at 6000 K (►) for a mixing ra¬ 
tio A^H 20 ’-^H 2 “32:32. 


From pairs of bound nuclei, we built larger chem¬ 
ical species if the bonds remained contiguous. Us¬ 
ing the first minimum in the radial distribution 
functions ( Soubiran &: Militzer ||2014 ), the follow¬ 
ing distance limits were derived: ^H-H = 1-0 A, 
Zq-h = 1-3 A and Iq-O = 1-9 A. 



Fig. 10.— Distribution of the hydrogen atoms 
among the different hydrogen bearing species at 
2000 K (□) and at 6000 K (►) for a mixing ratio 
-^H20:-^H2 =32:32. 

The Fig. [9| and [TO] show examples of the chemi¬ 
cal compositions of the mixture. We observed that 
H, H 2 , O, HO and H 2 O make up more than 99% 
of the observed species in the system. A few tran¬ 


sient bigger molecules were identified also. The 
temperature has a strong impact on the degree of 
dissociation as expected. While at 2000 K and 
below, the mixture is almost fully molecular, at 
6000 K we observe a system that is nearly fully 
dissociated. We thus expect a much higher electric 
conductivity at 6000 K because the dissociation is 
generally associated with ionization especially if 
hydrogen is present (Collins et aZ.||200T|. 


For both temperatures, the fraction of disso¬ 
ciated species appears to increase with pressure. 
Water molecules tend to split into hydroxide and 
hydrogen ions or even to fully dissociate, entering 
into the regime where no stable chemical bonds 
exist. 


4. Discussion 


The computed equation of state of the homoge¬ 
neous mixtures shows only small deviations from 
the ideal mixing approximation. We thus conclude 
that an ideal mixing law may provide quite reason¬ 
able estimates for the purpose of planetary mod¬ 
eling in general as long as no phase separation is 
expected and that very accurate EOSs for the pure 
systems are used. 


Moreover, we do not find a phase separation 
for the hydrogen-water mixtures we studied. As 
was observed for other mixtures at higher pressure 


(Wilson & Militzer|20I2a 

Gonzalez-Cataldo et al. 

2014 

1, the entropy is nearly ideal and contributes 


the most to the free energy of mixing, which sta¬ 
bilizes the homogeneous mixture. 


We have to stress here that our miscibility pre¬ 
dictions differ slightly from what has been in¬ 
ferred from high pressure experiments by |Bali et| 
al. ( 2013| ) (see Fig. for comparison). In these ex¬ 


periments, samples of different minerals were com¬ 
pressed up to a 2.5 GPa in presence of water. Due 
to chemical reactions, hydrogen was released. In 
some conditions a phase separation of water and 
hydrogen occurred, even for temperatures above 
1000 K. While these experiments may be relevant 
for the Earth’s interior, the conditions are quite 
different from the envelope of the giant planets 
where no minerals are present. We therefore sug¬ 
gest that diamond anvil cell experiments water- 
hydrogen mixtures in a mineral-free environment 
should be performed. 


With knowledge of the experimental work by 
































Bali et al. (2013), Aranovich (2013) constructed a 


van Laar mixing model using the molar volume of 
the pure species and a van Laar parameter fitted 
on earlier experiments at the kbar regime. This 
model predicts hydrogen-water mixture to phase 
separate at even higher temperatures that sug¬ 


gested by Bali et al. (2013). When we constructed 


a van Laar mixing model, we were not able to ob¬ 
tain a good fit to our ab initio Gibbs free energies, 
even though the molar volumes in our simulations 
of the pure species agree quite well with the values 
that Aranovich used. To match our ah initio data, 
the van Laar model would need to be extended by 
including nonideal mixing effects for the volume 
and the van Laar mixing parameter W would need 
to be made temperature and pressure dependent. 

The fact that our results indicate the water- 
hydrogen mix at 2-70 GPa has implications on our 
understanding of ice giant planets. The generic 
three layer model ( Nettelmann et aZ.|2013 ) with a 
hydrogen-rich outer envelope, a water-rich inter¬ 
mediate layer, and a rocky core may be a simplifi¬ 
cation. The boundary between the hydrogen and 
water layers is at about 10 GPa and 2000 K in 
Uranus and Neptune, which falls into the parame¬ 
ters regime that we explored with our simulations. 
Yet our results show that such a sharp boundary 
between the two layers is thermodynamically un¬ 
stable. Convection will efficiently mix the two lay¬ 
ers unless the system approaches a semi-convective 
state. 

Therefore we cannot rule out the three-layer 
picture completely because the time scale of the 
mixing has to be taken into account. In the core¬ 


accretion model (Pollack et al. 1996), solid plan- 


etesimals are accreted first, followed by a run-away 
gas accretion. The core is assumed to differenti¬ 
ate quickly letting the lighter components rise to 
the surface of the core. Based on our results, wa¬ 
ter would start to mix with the gas but the time 
scale is then a key parameter. If the mixing is only 
diffusive, the planet remains in a nearly differenti¬ 
ated state because the diffusive time scale of water 
in hydrogen can be as long as 10^^ years based on 
the diffusion coefficient estimates in ISoubiran fcl 


Militzer (2014) and the size of the envelope given 


by Nettelmann et al. (2013). On the other hand. 


if there is a vigorous and sustained convective ac¬ 
tivity then in a few convective time scale (on the 
order of a 100 years time scale (]Hubbard 1984, p 


294)) the planet should be homogenized, which 
would be incompatible with the multiple layers 
assumption. Nevertheless, [Leconte &: Ghabrier| 


(2012, 2013) showed that depending on the condi¬ 


tions, semi-convection may set in. If a gradient of 
composition exists in the interior of a planet, grav¬ 
ity may prevent an efficient convection of heavy 
elements. In this case, a series of alternating dif¬ 
fusive and convective layers are predicted to oc¬ 
cur. Although the long term evolution of such a 
semi-convective state is not fully understood it is 
a possible mechanism for maintaining fairly steep 
compositional gradients. The planet would then 
have a water-rich deep inner envelope and a grad¬ 
ual, semi-convective transition to a hydrogen-rich 
envelope. 

Finally, we want to discuss the importance of 
the water-hydrogen miscibility for exoplanet inte¬ 
rior modeling. Given a mass and radius obser¬ 
vation for a sub-Neptune planet we find the in¬ 
ferred hydrogen fraction depends significantly on 
whether water and hydrogen occur in mixed form. 
In Fig. [D we compare the interior properties of 
hypothetical planets composed of water and hy¬ 


drogen only. While Valencia et al. (2010) found 


very similar mass-radius relationships for homoge¬ 
neously mixed and differentiated, two-layer plan¬ 
ets composed of iron and silicates, here we find 
much larger deviations because hydrogen and wa¬ 
ter have very different compressibilities. For a 
planet of 10 Earth masses, we find the radius varies 
between 2.6 and 10.5 Earth radii for a pure water 
and pure hydrogen planets. Respectively, the cen¬ 
tral pressure varies between 5.8 and 0.07 Mbar. 

In Fig. El we also compare the hydrogen 
fraction for fully mixed and differentiated water- 
hydrogen planets for a given radius. If a planet 
with 10 Earth masses and 4 Earth radii were de¬ 
tected, one would infer a hydrogen fraction of only 
8% (0.8 Earth masses) if one assumed a differen¬ 
tiated interior. The central pressure would be 
5.6 Mbar, and the pressure at the water-hydrogen 
boundary would be 0.12 Mbar, which is far be¬ 
low the molecular-to-metallic transition pressure 
in pure hydrogen (Vorberger et al. 2007). On 


the other hand if one assumed a homogeneously 
mixed interior structure, the hydrogen mass frac¬ 
tion would increase by 25% (2.5 Earth masses) and 
the central pressure would decrease to 1.4 Mbar, 
which is sufficiently high for hydrogen molecules to 
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Fig. 11.— Predictions from mixed and differen¬ 
tiated interior models for H 2 -H 2 O planets with 10 
Earth masses. The upper panel shows the inferred 
H 2 mass fraction for a given radius. The lower 
panel displays the central pressure and, for differ¬ 
entiated planets, also the pressure at the boundary 
between the two layers. 


dissociate. This reaction would also increase the 
electrical conductivity of the mixture and further 
the generation of magnetic fields. 

The inferred hydrogen fraction increases signif¬ 
icantly if a mixed interior is assumed because hy¬ 
drogen gas is very compressible and then some hy¬ 
drogen fluid is exposed to much higher pressure. 
For differentiated planets of 10 Earth masses, we 
find hydrogen is never exposed to more than 0.31 
Mbar regardless of the planet’s radius. The cen¬ 
tral pressure was found to decrease with increasing 
hydrogen content because water is diluted with a 
material of smaller density. 

A difference of 0.8 and 2.5 Earth masses in hy¬ 
drogen contents also has implication for our un¬ 
derstanding of the environment where the planet 
formed. 

5. Conclusions 


the water-hydrogen mixtures behave close to an 
ideal mixture. We found that no phase separa¬ 
tion is expected for this parameter range. A ho¬ 
mogeneous mixture is always thermodynamically 
preferred for all concentrations. Our simulation 
results are in disagreement with water-hydrogen 


mixtures experiments in mineral cells (Bali et al. 


20131 in the lowest pressure-temperature regime 
explored by our calculations. We would suggest 
that high-pressure experiments of hydrogen-water 
mixture should be performed in a mineral-free en¬ 
vironment in order to better constrain the phase 
diagram in this regime. 

Since we predict no phase separation for wa¬ 
ter and hydrogen, this has consequences for the 
ice giant planets. If a planet already has a mixed 
water-hydrogen layer, the thermodynamic proper¬ 
ties, in particular the entropy of mixing, prevent 
any differentiation from occurring. 


If a planet has two separate water and hydro¬ 
gen layers, this implies two things. First, the icy 
materials must have been delivered early on when 
very little gas was present or the icy planetesimals 
must have been sufficiently large so that they pen¬ 
etrated through the existing atmosphere. Both 
possibilities would be consistent with the core ac¬ 
cretion model. 


Furthermore such a planet could not be fully 
convective, otherwise water and hydrogen would 
mix assuming the pressure at the interface is in 
the 2-70 GPa range. It is possible, however, for 
the planet to remain predominantly differentiated 
because the interior may assume a semi-convective 
state, in which water and hydrogen would not mix 
efficiently. Because the long-term dynamics of this 
semi-convection is not sufficiently well understood, 
it is difficult to model the evolution of planetary 
interiors on the gigayear time scale when compo¬ 
sitional gradients are present. 

We also showed that, unlike in the case of rock- 
iron mixtures, the mixing of water and hydrogen 
can drastically increase the estimates of the hy¬ 
drogen content of sub-Neptune exoplanets. This 
effect has to be taken into account when their in¬ 
terior structure and evolution are modeled. 


Using DFT molecular dynamics simulations, we 
showed that on the range of pressure from 2 to 
70 GPa and temperature from 1000 to 6000 K, 


10 













Acknowledgments 


This work was supported by NASA and NSF. 

Computers at NAS and NCOS were used. 

REFERENCES 

Aranovich, L. Ya. 2013, Petrology, 21, 539. 

Baldereschi, A. 1973, PhRvB, 7, 5212. 

Bali, E., Audetat, A. & Keppler,H. 2013, Natur, 
495, 220. 

Blochl, P. E. 1994, PhRvB, 50, 17953. 

Bodenheimer, P. & Lissauer, J. J. 2014, ApJ, 791, 
103. 

Caillabet, L., Mazevet, S. & Loubeyre, P. 2011, 
PhRvB, 83, 094101. 

Cavazzoni, C., Chiarotti, G. L., Scandolo, S. & 
Tosatti, E. 1999, Sci, 283, 44. 

Collins, L. A., Bickham, S. R., Kress, J. D., 
Mazevet, S., Lenosky, T. J., Trouiller, N. J. & 
Windl, W. 2001, PhRvB, 63, 184101. 

Connerney, J. E. P., Acuna, M. H. & Ness, N. F. 
1987, JGR, 92, 15329. 

Connerney, J. E. P., Acuna, M. H. & Ness, N. F. 
1991, JGR, 96, 19023. 

Dion, M., Rydberg, H., Schroder, E., Langreth, 
D. G. & Lundqvist, B. 1. 2004, PhRvL, 92, 
246401. 

Encrenaz, T. 2003, P&SS, 51, 89 

Exoplanet Archive, URL: 

ht tp: / / exoplanet archive. ipac. caltech.edu, 

Fegley, B. & Prinn, R. G. 1986, ApJ, 307, 852. 

Fortney, J. J. & Hubbard, W. B. 2004, ApJ, 608, 
1039. 

Fortney, J. J., Mordasini, C., Nettelmann, N., 
Kempton, E. M.-R., Greene, T. P. & Zahnle, 
K. 2013, ApJ, 775, 80. 

French, M., Mattsson, T. R., Nettelmann, N. & 
Redmer, R. 2009, PhRvB, 79, 054107. 


French, M. & Redmer, R. 2009, JPGM, 21, 
375101. 

Gonzalez-Gataldo, F., Wilson, H. F. & Militzer, 

B. 2014, ApJ, 787, 79. 

Guillot, T. 1999, Sci, 286, 72. 

Guillot, T. 2005, AREPS, 33, 493. 

Helled, R., Anderson, J. D., Podolack, M. & Schu¬ 
bert, G. 2011, ApJ, 726, 15. 

Helled, R. & Bodenheimer,P. 2014, ApJ, 789, 69. 

Hohenberg, P. & Kohn, W. 1964, PhRv, 136, 864. 

Hubbard, W. B. 1968, ApJ, 152, 745. 

Hubbard, W. B., Planetary Interiors, Van 
Norstrand Reinhold, 1968. 

Hubbard, W. B., Podolack, M. & Stevenson, D. J. 
1995. The interior of Neptune. In: Gruikshank, 
D. P. (Ed.), Neptune and Triton. University of 
Arizona Press, Tucson, pp. 109-138. 

Hubbard, W. B. 1999. Interiors of the giant plan¬ 
ets. In: Beatty, J. K., Petersen, G. C., Ghaikin, 
A. (Eds.), The New Solar System. Sky Publish¬ 
ing Go., Gambridge, MA, pp. 193-200. 

Jacobson, R. A. 2007, BAAS, 39, 453. 

Jacobson, R. A. 2009, AJ, 137, 4322. 

Klimes, J., Bowler, D. R. & Michaelides, A. 2010, 
PhRvB, 83, 195131. 

Knudson, M. D., Hanson, D. L., Bailey, J. E., Hall, 

C. A., Asay, J. R. & Deeney, G. 2004, PhRvB, 
69, 144209. 

Knudson, M. D., Desjarlais, M. P., Lemke, R. W., 
Mattsson, T. R., French, M., Nettelmann, N. & 
Redmer, R. 2012, PhRvL, 108, 091102. 

Kohn, W. & Sham, L. J. 1965, PhRv, 140, 1133. 

Kresse, G. & Furthmiiller, J. 1996, PhRvB, 54, 
11169. 

Leconte, J. & Ghabrier, G.2012, A&A 540, A20. 

Leconte, J. & Ghabrier, G.2013, NatGe 6, 347. 

Lenosky, T. J., Kress, J. D. & Gollins, L. A. 1997, 
PhRvB, 56, 5164. 


11 



Loubeyre, P., Brygoo, S., Eggert, J., Celliers, 
P. M., Spaulding, D. K., Rygg, J. R., Boehly, 
T. R., Collins, G. W. & Jeanloz, R. 2012, 
PhRvB, 86, 144115. 

McMahon, J. M., Morales, M. A., Pierleoni, C. & 
Ceperley, D. M. 2012, RvMP, 84, 1607. 

Mermin, N. D. 1965, PhRv, 137, A1441. 

Militzer, B., Hubbard, W. B., Vorberger, J., Tam- 
blyn, I. & Bonev, S. A. 2008, ApJ, 688, L45. 

Militzer, B. 2013, PhRvB, 87, 014202. 

Monkhorst, H. J. & Pack, J. D. 1976, PhRvB, 13, 
5188. 

Morales, M. A., Pierleoni, C., Schwegler, E. & 
Ceperley, D. M. 2009, PNAS, 106, 1324. 

Ness, N. E., Acuna, M. H., Behannon, K. W., 
Burlaga, L. E., Connerney, J. E. P., Pepping, 
R. P. & Neubauer, F. M. 1986, Sci, 233, 85. 

Ness, N. F., Acuna, M. H., Burlaga, L. F., Con¬ 
nerney, J. E. P., Pepping, R. P. & Neubauer, 
F. M. 1986, Sci, 246, 1473. 

Nettelmann, N., Helled, R., Fortney, J. &Redmer, 
R. 2013, P&SS, 77, 143 . 

Nose, S. 1984, JChPh, 81, 511. 

Nose, S. 1991, PThPS, 103, 1. 

Pan, D., Spanu, P., Harrison, B., Sverjensky, D. A. 
& Gain, G. 2013, PNAS, 110, 6646. 

Pearl, J. C., Conrath, B. J., Hanel, R. A., Pir- 
raglia, J. A. & Coustenis, A. 1990, Icar, 84, 12. 

Pearl, J. C. & Conrath, B. J. 1991, JGRS, 96, 
18921. 

Perdew, J. P., Burke, K. & Ernzerhof, M. 1996, 
PhRvP, 77, 3865 . 

Petigura, E. A., Howard, A. W. & Marcy, G. W. 
2013, PNAS, no, 19273 

Podolak, M., Reynolds, R. T. & Young, R. 1990, 
GeoRP, 17, 1737. 

Pollack, J. B., Hubickyj, O., Bodenheimer, P., 
Pissauer, J. J., Podolak, M. & Greenzweig, Y. 
1996, Icar, 142, 62. 


Santra, B., Klimes, J., Tkatchenko, A., Alfe, D., 
Slater, B., Michaelides, A., Car, R. & Scheffler, 
M. 2013, JChPh, 139, 154702. 

Seager, S., Kuchner, M., Hier-Majumder, C. A. & 
Militzer, B. 2007, ApJ, 669, 1279. 

Seward, T. M. & Franck, E. Pl. 1981, Ber. Bunsen 
Phys. Chem., 85, 2. 

Soubiran, F., Mazevet, S., Winisdoerffer, C. & 
Chabrier, G. 2013, PhRvB, 87, 165114. 

Soubiran, F. & Militzer, B. 2014, HEDP. 

Stevenson, D. J. & Salpete, E. E. 1977, ApJS, 35, 
239. 

Stevenson, D. J. 1985, Icar, 62, 4 

Tamblyn, P, Bonev, S. A. 2010, PhRvP, 104, 
065702. 

Valencia, D., Ikoma, M., Guillot, T. & Nettel¬ 
mann, N. 2010, A&A, 516, A20. 

Vorberger, J., Tamblyn, P, Militzer, B. & Bonev, 
S. A. 2007, PhRvB, 75, 024206. 

Wahl, S. M., Wilson, H. F. & Militzer, B. 2013, 
ApJ, 773, 95. 

Wahl, S. M., Wilson, H. F. & Militzer, B. 2015, 
E&PSP, 410, 25. 

Weiss, P. M. & Marcy, G. W. 2014, ApJ, 783, P6. 

de Wijs, G. A., Kresse, G. & Gillan, M. J. 1998, 
PhRvB, 57, 8223 

Wiktorowicz, S. J. & Ingersoll, A. P. 2007, Icar, 
186, 436 

Wilson, H. F. & Militzer, B. 2010, PhRvP, 104, 
121101 . 

Wilson, H. F. & Militzer, B. 2012, ApJ, 745, 54. 

Wilson, H. F. & Militzer, B. 2012, PhRvP, 108, 
111101 . 

Wilson, H. F., Wong, M. P. & Militzer, B. 2013, 
PhRvP, no, 151102. 

Wong, M. H., Mahaffy, P. R., Atreya, S. K., Nie¬ 
mann, H. B. & Owen, T. C. 2004, Icar, 171, 
153. 


This 2-column preprint was prepared with the AAS 
macros v5.2. 


12 



